Gauge Drivers for the Generalized Harmonic Einstein Equations 
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The generalized harmonic representation of Einstein's equation is manifestly hyperbolic for a 
large class of gauge conditions. Unfortunately most of the useful gauges developed over the past 
several decades by the numerical relativity community are incompatible with the hyperbolicity of 
the equations in this form. This paper presents a new method of imposing gauge conditions that 
preserves hyperbolicity for a much wider class of conditions, including as special cases many of 
the standard ones used in numerical relativity: e.g., -freezing, F-freezing, Bona-Masso slicing, 
conformal F-drivers, etc. Analytical and numerical results are presented which test the stability and 
the effectiveness of this new gauge driver evolution system. 



I. INTRODUCTION 

The gauge (or coordinate) degrees of freedom in the 
generahzed harmonic (GH) form of the Einstein equa- 
tions are determined by specifying the gauge source func- 
tions H"^ . These functions are defined as the results of 
the covariant scalar-wave operator acting on each of the 
spacetime coordinates x": 



(1) 



The GH form of the Einstein equations can be repre- 
sented (somewhat abstractly) as 

V'"'9c5<Jl/.ab -I- daHt, + dbHa = Qab{H, V, ^V"), (2) 

where ipab is the spacetime metric, Ha = ipabH^ i and 
Qab represents lower order terms that depend on Ha , the 
metric, and its first derivatives. 

The GH form of the Einstein equation is manifestly 
hyperbolic whenever Ha is specified as an explicit func- 
tion of the coordinates and the metric: Ha = Ha{x,ip). 
In this case the terms daHi, that appear in Eq. ^ con- 
tain at most first derivatives of the metric. The Einstein 
equations become, therefore, a set of second-order wave 
equations for each component of the spacetime metric: 



Ip" dcddlpab = Qabix, dip). 



(3) 



Thus the Einstein equations are manifestly hyperbolic for 
any Ha = Haix,i>). 

Most of the useful gauge conditions developed by the 
numerical relativity community over the past several 
decades can not, unfortunately, be expressed in the sim- 
ple form Ha = Ha(x, ip) (unless the full spacetime metric 
'>Pab = V'ab[a;] IS known a priori). Many of these condi- 
tions (e.g., maximal slicing or F-drivers) would require 
gauge source functions that depend on the spacetime 
metric and its first derivatives: Ha — Ha{x,^,d'ip). In 
this case the terms daHb in Eq. ([2]) depend on the sec- 
ond derivatives of the metric, if^ab, and this (generically) 
destroys the hyperbolicity of the system. 



Pretorius [l|, 0, 0] proposed a way to expand signifi- 
cantly the class of allowed gauge conditions, by elevating 
Ha to the status of an independent dynamical field. A 
separate gauge driver equation is introduced to evolve 
Ha, for example. 



WcHa^Qa{x,H, dH,^,d^P), 



(4) 



where V^VciJa denotes the wave operator^ acting on Ha- 
This gauge driver equation is solved together with the 
GH Einstein equations to determine ipab and Ha simul- 
taneously. In the combined evolution system, consisting 
of Eqs. (H]) and (g]), the daHb terms in Eq. ^ are now 
lower-derivative terms that do not affect the hyperbolic- 
ity of the system. Thus the combined GH Einstein plus 
gauge driver system is manifestly hyperbolic so long as 
Qa on the right in Eq. (g]) depends only on the fields and 
their first derivatives: Qa = Qa{x, H,dH,^,dip). Each 
of the solutions, Ha — Ha{x), to these gauge driver equa- 
tions is a gauge condition. So the gauge driver system 
provides a way for Ha to be determined by the metric 
and its derivatives in a flexible way without destroying 
the hyperbolicity of the GH Einstein equations. 

Pretorius used a particular gauge driver equation of 
this form to determine Ht in his ground breaking binary 
black-hole simulations: 

W^W.Ht ^ Ci{l~N)N-PH2{dtHt-N''dkHt)N-\ (5) 

where in this case V^Vc is the covariant scalar-wave op- 
erator, N is the lapse, TV*"' is the shift, and ^i, ^2 and p 
are constants. For suitable choices of these parameters, 
Pretorius found this system to be quite effective in pre- 
venting the lapse from "collapsing" toward zero as the 
system evolves. Solutions to this gauge driver do not 
correspond to any of the traditional gauge conditions of 
numerical relativity as far as we know. 



^ We define exactly what we mean by this wave operator in 
SecHlBl 
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In this paper we introduce a new class of gauge driver 
equations that are general enough to provide implemen- 
tations of (almost) all of the standard gauge conditions 
used by the numerical relativity community. This is 
done by choosing an appropriate "source" term Qa for 
the right side of Eq. The idea is to choose Qa so 
that solutions Ha evolve quickly toward a target gauge 
source function Fa ■ This Fa is chosen so that strict equal- 
ity Ha = Fa corresponds exactly to the gauge condition 
of interest to us. We limit Fa only by assuming that 
it depends on the spacetime metric and its first (but 
not second) derivatives: Fa — Fa{x,tp,dtjj)- These new 
gauge driver equations are introduced in Sec. |TT1 and we 
show there that the combined GH Einstein plus gauge 
driver system is symmetric hyperbolic for any target 
gauge source function of this allowed form. In Sec. IIIII 
we present the target gauge functions Fa correspond- 
ing to (many of) the gauge conditions commonly used 
by the numerical relativity community, including max- 
imal slicing, iiT-freezing, Bona-Masso slicing, conformal 
F-freezing, and conformal F-drivers. In Sec. IIVI we use 
analytical methods to analyze the solutions of the new 
gauge driver system. We show in particular that Ha 
approaches any (time independent) Fa exponentially for 
evolutions of the gauge driver equations on flat space. 
We also demonstrate the effectiveness of the coupled GH 
Einstein and gauge driver system for the case of small 
perturbations of flat space using Bona-Masso slicing and 
one of the conformal F-driver conditions. In Sec. IVl we 
show the effectiveness and stability of our implementa- 
tion of a particular choice of Bona-Masso slicing and con- 
formal F-driver condition using numerical solutions of the 
full non-linear equations for perturbed single black-hole 
spacetimes. We summarize and discuss these various re- 
sults in Sec. EH 



II. GAUGE DRIVER EQUATIONS 

We begin this section by deriving a system of gauge 
driver equations in Sec. Ill Al and then constructing a 
general first-order representation of these equations in 
Sec. IIIBI We derive the characteristic fields for this sys- 
tem and their associated speeds in Sec. Ill C[ and show 
that the coupled gauge-driver and GH Einstein system 
is symmetric hyperbolic. We analyze the constraints in 
Sec. Ill D1 and derive constraint preserving boundary con- 
ditions for the gauge-driver fields in Sec. Ill El 



A. Motivation 

In this section we provide some motivation for our 
choice of gauge driver equation. We consider first the 
case of the gauge driver V^\7cHa — Qa acting on a fixed 
flat-space background. The idea is to choose Qa so that 
the solutions, Ha-, to this equation quickly approach the 
desired target gauge source function Fa- If Fa were con- 



stant in space and time, there would be a fairly obvious 
and simple choice: 

WcHa=Qa= f^\Ha-Fa)+2ndtHa, (6) 

where ^ is a freely specifiable constant. If Ha like Fa 
were independent of spatial position, the gauge driver 
equation would be equivalent in this case to the ordinary 
differential equation, 

dfiHa - Fa) + 2^ldt{Ha ~ Fa) + fl^iHa - Fa) - 0, (7) 

whose solution has the form: Ha{t) — Fa + [Ha{0) — 
^a]e~'^*. A similar argument applied to the spatial 
Fourier transform of Eq. ([6]) shows that spatially inho- 
mogeneous Ha also approach Fa exponentially in time. 
In this special case (i.e., spatially homogeneous and time 
independent Fa) the simple gauge driver has the desired 
behavior: all the solutions Ha approach the target gauge 
source Fa exponentially on the adjustable time scale 

This simple gauge driver, Eq. ([5]), fails unfortunately 
even in flat space if Fa is a generic function of position. 
An easy way to see this is to assume that all the solu- 
tions Ha do approach Fa asymptotically as t — > oo. Since 
Fa and consequently Ha are independent of time in this 
limit, Eq. © reduces to V^VkHa = 0, where V'^VuHa 
represents the spatial Laplacian of Ha- But this is im- 
possible because ^^^kFa ^ for generic Fa- So (not 
surprisingly) the simple gauge driver fails in general. 

This gauge driver can be modified in a fairly straight- 
forward way, however, that corrects this problem. Define 
an auxiliary dynamical field 9a'- 

dtOa+'nOa^V'^VkHa- (8) 

This equation can be integrated analytically to obtain an 
equivalent integral representation of 9a'- 

Oa = ea(0)e-"* + /e-''^*-*') V'^Vfei/,(i') dt' - (9) 

Thus 6a represents an exponentially weighted (in favor 
of times near t) time average of the past evolution of the 
term on the right side of Eq. ([5]) . We can use this 9a to 
construct an improved gauge driver: 

W,Ha =Qa= ^i^{Ha - Fa) + 2^ldtHa + ^Oa- (10) 

If a solution to the improved gauge driver approaches 
a time independent state, then Eq. ([8]) implies that 
i]9a = V^VkHa- Equation pH)) reduces in this case to 
— fJ,^{Ha — Fa)- So the addition of the time averaging 
field 9a forces Ha to approach Fa in any time indepen- 
dent state, even in the case of inhomogenous Fa- The 
remainder of this paper is devoted to the analysis of this 
improved gauge driver, Eq. pop . suitably generalized for 
use in an arbitrary spacetime. 
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B. First Order Form 

We find it very useful to consider the first-order repre- 
sentations of evolution systems (such as our gauge driver) 
for a variety of reasons: from basic mathematical issues 
(such as the formulation of appropriate boundary con- 
ditions), to more practical code stability issues. This 
section develops a first-order representation of the gauge 
driver system, suitably generalized for use in an arbitrary 
curved spacetime. 

The gauge driver system described above evolves Ha 
through a wave equation of the form V^X/cHa = Qa- 
In a general curved spacetime, we assume that V^VciJa 
represents the covariant wave operator that treats Ha 
as a co-vector. This choice needs a bit of clarification, 
since the gauge source function Ha is not actually a co- 
vector. One way of giving meaning to this equation is 
to use the gauge driver system to determine a new field 
Ha that does transform as a co- vector: V^V cHa — Qa- 
Then fix Ha by setting Ha = Ha in some particular co- 
ordinate frame. This construction is not covariant, but 
fixing coordinate conditions can never be completely co- 
variant. An equivalent way to do this is to write out 
and impose the gauge driver equation, V^Vci?a — Qa, 
only in the special coordinate frame in which Ha = Ha- 
We adopt this second approach since it simplifies the no- 
tation somewhat. So throughout this paper the gauge 
driver equation, VVc^^a — Qa, will only be imposed in 
some particular coordinate frame that we must specify. 
In our code we use a global Cartesian coordinate system, 
and we will always impose the gauge driver equation in 
that frame. 

Before we discuss the first-order form of the gauge 
driver equations, we also need to examine the somewhat 
pathological covariant vector wave operator in more de- 
tail. This operator acting on Ha (assumed here to be a 
co-vector as discussed above) can be written out more 
explicitly in the form: 

v^VeiJa = tp''^dbdcHa-r'dbHa-2i;''-ri^dbHd 

+ {Ra'' - daT'')Hb, (11) 

where F'^bc is the Christoffel connection, = -^^^F^tc, 
and Ra'^ is the associated Ricci curvature. This wave op- 
erator is well behaved on a fixed background spacetime. 
However the H^daT'' term includes second derivatives of 
the metric that would interfere with hyperbolicity, if it 
were coupled in a non-trivial way to the full Einstein 
equations. Fortunately this problem has a simple solu- 
tion. Since we use the GH form of the Einstein equations, 
this term can be transformed into the more benign form, 
—HbdaH'' (or if a more linear looking form is preferred 
TbdaH''), using the gauge constraint H"^ — —T°- [4j. We 
regard the Ricci tensor Ra'' as being determined by the 
matter sources via the Einstein equations; in particular, 
it does not contain any second derivatives of the met- 
ric. For notational convenience we introduce the quantity 



Wa{H), 

WaiH) = 2ij'''T^acdbHd-{daH'' + Ra'')Hb, (12) 

that represents the parts of the vector wave operator that 
are not present in the scalar wave operator. Our repre- 
sentation of the covariant vector wave operator is there- 
fore given by, 

W,Ha = ^'''dbd,Ha~T''dbHa-Wa{H). (13) 

To represent this equation in first-order form, wc in- 
troduce the usual additional first-order fields 11^ and 
representing (up to the addition of constraints) the ap- 
propriate time and space derivatives of Ha respectively: 

nf = -t'dbHa, (14) 
'^fa = d.^Ha- (15) 

Here (and throughout this paper) t° is the future directed 
unit normal to the t = constant hypersurfaces; Latin 
indices a through h are spacetime indices and run from 
to 3; and Latin indices i through n are spatial indices 
and run from 1 to 3. We also define the spatial metric 
on the t = constant hypersurfaces, 

Qab^ i'ab + tah- (16) 

The covariant wave operator, VVcHa, can then be ex- 
pressed in terms of these first-order field variables: 

V^V.Ha - t^dcUf +g'^d,^fa-tbT'n^ -T'P^a 

-finft^i^Hfe^ + g'^^^t^n,, 

-Wa{H), (17) 
where Wa{H) can be written as 

Wa{H) = (talibc+ ga'^^bc)t^r''{^d ~t"HdH,) 
Hta^^b + 9a'^J^b)^*'•'9'\^kc + HuH,) 

-g'H''n,a^% + g'^i;''-'i>,ab^f^ 
-5^mf (H„ -t- t'><^,ab) - .g'^g'^'$,fca$f^ 

+ {talif + ga^'^fb)^'' - Ra^'Hb- (18) 

We note that leaving out the Wa{H) terms is equivalent 
to applying the covariant scalar wave operator to each 
component of Ha in our special coordinate frame. We 
also note that the remaining F** terms that appear in the 
above equations are to be thought of as functions of the 
first-order GH fields: 

F" = ^"^i'^H.rf -I- 5^^ V^'^<i>yc - \i^'\t''n,a + S^'^^cd). 

(19) 

The representation of wave equations of this type in 
first-order form is well understood, see e.g., Refs. 0, [H]; 
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the result for our gauge driver equation is 

dtHa - (1 + lf )N^dkHa - - A^nf - 7f iV"<i>fa, 

(20) 

+ Qa+NWa, (21) 

-li^d.N + <i>f,9,7V'= - 72"^$^I- (22) 

The quantities iV, iV*"', and gij that appear in these equa- 
tions are the lapse, shift and spatial metric, defined by 
the usual three-plus-one representation of the spacetime 
metric: 

ds^ = tpabdx"-dx'' 

= -N^dt^ + g,j{dx' + N'dt){dx^ + N^dt). (23) 

The auxiliary quantities K and J* in Eq. (|2ip depend on 
the background spacetime geometry and can be written 
in terms of the first-order GH Einstein variables: 

K = iff'^y (24) 

J' = (.g^ V' - ^5^^/') <i>,M + ^5^^i"i%a6. (25) 

The constants 7f and ^2 introduced (in analogy with 
the first-order GH Einstein system [4]) to allow us to 
control the growth of constraint violations, and to allow 
us to adjust one of the characteristic speeds of the system. 

The quantity Qa in Eq. (|2ip is defined by a natural 
generalization of Eq. (fTO|) : 

Qa = ^il{l-il)NiHa-Fa) 

-2/i2(l-6)^nf + 7710,. (26) 

The differences between this expression and Eq. pO)) are 
an overall factor of the lapse N (to convert from coordi- 
nate time to proper time), the replacement of dtHa by 
(the first order field representing — f^dcHa), the in- 
troduction of independent damping parameters fii and 
fi2, and the introduction of new constant parameters 
and ^2- The purpose of these latter parameters, 
^1 and ^2, is to move the damping terms (or fractions 
thereof) into the source for the time-averaging field da (cf. 
Eq. (ITT)) below) , thus effectively replacing these terms by 
their time averages. We assume as before that Fa is a 
given function of the four-metric and its first derivatives: 

Fa=Faix,^,d^P). 

The evolution equation for 9a is chosen, in analogy 
with Eq. ([5]), to include as its source all the terms in 
Eq. (j21[) that do not vanish automatically in a time in- 
dependent state: 

dtOa + VlOa = 2M2K3A^nf + (1 - + -f^)N^dkHa] 

-2^*2(1 - 6)7f + 7f 72"^^'*fa 

-ivisrnf - Nr^fa - NWa 

+fil^iN{Ha - Fa) - 2/^26A^nf . (27) 



The ^3 parameter is introduced to add a multiple of dtFla 
to the source of the time averaging field. We use Eq. (^0]) 
to re-express this dtHa as the terms proportional to ^3 
that appear on the right side of Eq. (P7)) . Assuming the 
system approaches a state in which Ha becomes time in- 
dependent, then da exponentially approaches the time in- 
dependent limit of the terms on the right side of Eq. (j27|) . 
These terms were chosen so that Eq. ([?T|) then implies 
that Ha — > Fa in this limit. Our choices for the param- 
eters, /^i, /i2, ?7i, Cii and ^3 that appear in Eqs. ([26|l 
and (|27p will be guided by the stability analysis that we 
perform in Sec. IIVI 

C. Characteristic Fields 

The gauge driver evolution Eqs. (I20l)-(l22l) and ((27|) 
comprise a first-order evolution system of the form, 

dtu'^ + A^'^ndkvp = B", (28) 

for the fields = {i?a, n^, Qa] (treating the space- 
time metric for the moment as a fixed background field). 
The characteristic fields of such an evolution system are 
important for a number of reasons, including the for- 
mulation of outer boundary conditions and exchanging 
information across internal boundaries of the computa- 
tional domain. The characteristic fields (in the direction 
of a unit spacelike covector n^) are defined as the pro- 
jections of the fields onto the left eigenvectors of the 
characteristic matrix UkJ^"" p. For the gauge driver sys- 
tem, these characteristic fields are 

;7f± = nf ±n'$^-7fiJ„ (29) 

= Ha. (30) 

Zfa^ = P^^l. (31) 

^ -2^2(1 -6)i?a, (32) 

where = gij — ninj. 

The eigenvalues associated with the characteristic 
fields are called the characteristic speeds of the system. 
For the gauge driver system, the characteristic fields Ua^ 
have speeds ±N — riiN'^, Z^^ has speed — (H- 7f )niiV*, 
Z^^^ has speed —riiN^, and Za^ has speed zero. 

The inverse transformation between dynamical fields 
and characteristic fields for our gauge driver system is 



Ha 


ryHl 

— : 




(33) 


n" 






(34) 




- \iu»* - 


-U^-)n., + Zg\ 


(35) 


Oa 




iUa^^ + Ua"-) 






+2Ai2(l- 


-^s)Zr~7iZ^\ 


(36) 



The existence of this inverse transformation shows that 
there is a one-to-one correspondence between the dynam- 
ical fields and the characteristic fields. This implies that 
the gauge driver system is strongly hyperbolic. 
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A quasi-linear evolution system, Eq. ((28)) . is symmetric 
hyperbolic (a stronger condition than strong hyperbolic- 
ity) if there exists a positive definite metric Sa/3 (called a 
symmetrizer) on the space of fields, such that Sa^A'^'^fS — 
Si^-yA'^'^a- The gauge driver system, Eqs. (I20l)-(l22l) and 
P7|) . does have such a symmetrizer: 



a 

+g'H^fj^f^ + (dnf - il^dHaf], (37) 

where are arbitrary (non- vanishing) constants. The 
gauge driver system is therefore symmetric hyperbolic. 

Up to this point the discussion has focused on the prop- 
erties of the gauge driver system, Eqs. ([20|) ~([22 |l and (pTjl . 
with the spacetime metric considered as a fixed back- 
ground field. Our real interest of course is the case where 
the gauge driver system is coupled to the GH Einstein 
system, Eq. ^ . Thus we need to consider the properties 
of the combined evolution system having as dynamical 
fields the gauge driver fields plus the GH Einstein system 

fields: = {i^a,^^,$,^,0„Vab,^ab,$^ab}. The ficlds 

Hab and <^iab represent the first derivatives of the space- 
time metric tpab, as defined for example in Ref. We 
need to analyze the properties of the characteristic ma- 
trix A^°' p of this combined system to determine whether 
the full coupled system is hyperbolic. 

We have shown above that the gauge driver system, 
Eqs. (HOll-jll]) and ((271), is symmetric hyperbohc if the 
spacetime metric is considered as a background field. 
Similarly, the first-order representation of the GH Ein- 
stein system is symmetric hyperbolic if the gauge 
source function Ha is considered as a background field. 
The characteristic matrix A^"' p of the combined system 
is block diagonal, except for any cross terms that might 
arise if derivatives of the gauge driver fields appear in evo- 
lution equations for the GH fields or vice versa. The only 
potential cross terms are as follows: The term 2d(^aHh) 
occurs in the GH Einstein equations ([2]), the quantities K 
and Ji in Eq. (PT|) depend on derivatives of the spacetime 
metric, and the target gauge source function Fa appear- 
ing in Eq. (|2ip may include derivatives of the spacetime 
metric. 

However, di^aHh) can be rewritten in terms of the first- 
order gauge driver variables as 



d(aHb) 



(a9b) 



(a^b)- 



(38) 



Likewise, K and Ji can be expressed as algebraic func- 
tions of the first-order GH fields ipab, ^ab, and ^iab, cf. 
Eqs. ([M]) and (P5|) . Finally, the target gauge source func- 
tion Fa is assumed to be a function of the metric and its 
first derivatives, so it also can be written as an algebraic 
function of the first-order fields: Fa — Fa{x,il;,Il,^). 
Thus all of these potential cross terms can be written 
as algebraic functions of the dynamical fields and do not 
contribute to the characteristic matrix A'^"^ at all. 



The characteristic matrix of the combined evolution 
system is therefore block diagonal. It follows that the 
characteristic fields of the combined system are just the 
collection of unmodified characteristic fields from the sep- 
arate systems. Similarly the matrix Saf3 needed to sym- 
metrize the full system is just the matrix whose diagonal 
blocks are the symmetrizers of the individual systems. 
It follows trivially that the combined GH Einstein and 
gauge driver system is both strongly and symmetric hy- 
perbolic. 



D. Constraints 

The basic gauge driver evolution system, Eq. (jH), has 
no fundamental constraints. However by transforming 
the system to first-order form, Eqs. we intro- 

duce a set of new constraints: 



Cii — diHa - <i>^. 



C 



H 



(39) 
(40) 



These constraints vanish, = Cf^„ = 0, if and only 
if a solution to the first-order system also represents a 
solution to the original second-order equation. 

These constraints are determined by the values of the 
dynamical fields Ha and therefore their time evolu- 
tion is determined by the gauge driver evolution system. 
It is straightforward to show that these constraints sat- 
isfy the evolution equations 



(1 + 7r)N^dkC = -(1 + 7f )9.A^'Cf„ 

-72"^C.^-7fA^'Ql, (41) 

-272«%ArCjf„ (42) 



as a consequence of Eqs. ([201) HEU)- 

The characteristic matrix of this constraint evolution 
system is diagonal, so the constraints are themselves 
characteristic fields of this system. The constraint C^^ 
propagates at the speed —{l + -f")nkN'^, while Cfja prop- 
agates at the speed —UkN^. This constraint evolution 
system is strongly (and also symmetric) hyperbolic. 

The constraint evolution system is also homogeneous 
in the constraints, i.e., the right sides of Eqs. ([H]) and 
(|42)) are proportional to the constraints. This implies, 
for example, that these constraints will remain satisfied 
within the domain of dependence of the subset of the 
initial surface on which they are satisfied. 



E. Boundary Conditions 

Boundary conditions are needed for any of the charac- 
teristic fields having incoming (i.e., negative) character- 
istic speeds on the boundary. Some of these boundary 
conditions can be determined by the need to prevent the 
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influx of constraint violations, while others can be cho- 
sen to control the particular gauge condition being im- 
posed at the boundary. In analogy with the scalar field 
system [B*], the needed constraint preserving boundary 
conditions for this system are: 

dtZ^' = DtZ^^ - (1 + 7f )nfc7V'=n'C,^, (43) 
dtZg-" = D^ZS^^-n^N^P^^u'Cfi,, (44) 

where dtZ^^ = dtHa and dtZf^"^ = Pi'^dt^j^^ represent 
the constraint field projections of the time derivatives of 
the dynamical fields, while DtZ^^ and DfZ^^ represent 
the constraint field projections of the right sides of the 
evolution equations for these fields. 

The characteristic fields U^^ need boundary condi- 
tions whenever the corresponding speeds v± — ±N — 
n}^N^ are negative. Since V- < typically the U^^ 
mode is the one needing a boundary condition. The 
boundary condition on this field controls the incom- 
ing part of the gauge condition being imposed on the 
boundary. We often use a "freezing" boundary condi- 
tion, dtU^^ = 0, or the boundary condition, dtU^" = 
— 7^9tZ^^. Ideally the boundary condition on the char- 
acteristic field should be determined by the gauge 
condition that the driver equation is trying to enforce, 
however at present we do not know how to do this. 



III. SPECIFIC GAUGE CONDITIONS 

The gauge driver equations presented in Sec. |TT] were 
designed to evolve the gauge source function. Ha, toward 
a target function Fa — Fa{x, ip, d^). The question of how 
well these equations accomplish this will be explored in 
Sees. HV] and |Vl Here we focus on the issue of construct- 
ing target functions Fa for particular gauge conditions 
used in numerical relativity. 

Most of the gauge choices used by the numerical rel- 
ativity community, including all the examples below, 
are expressed as conditions on the spacetime metric 
and its first (space and time) derivatives; so abstractly, 
all such gauge conditions can be written in the form 
Ga{x, ip, dtp) = 0. Whenever the GH Einstein constraints 
are satisfied, it follows from Eq. ^ that Ha — —Ta = 
—Tabc'4''"^j where Tate is the four-dimensional Christoffel 
symbol. An appropriate target gauge source function Fa 
is therefore given by 

Fa = -Ta-pGa, (45) 

where p is an arbitrary (non- vanishing) constant. When 
the constraints are satisfied, this equation implies that 
Ha — Fa — pGa- So if the gauge driver system succeeds 
in driving Ha — Fa ^ 0, it follows that Ga —> as well 
for any p ^ Q. This Fa has the general form assumed 
in the discussions of Sec. HH i^a = Fa{x , tp , dip) , when- 
ever Ga has the form Ga = Ga{x,tp,dip). Therefore the 
gauge driver system with this target Fa should enforce 



the desired gauge condition Ga = asymptotically as the 
system evolves. 

The numerical relativity community traditionally sep- 
arates gauge conditions into those that determine the 
lapse N (often called slicing conditions) and those that 
determine the spatial coordinates through the shift N'^. 
Expressing Ta in terms of the three-plus-one representa- 
tion of the spacetime metric, Eq. reveals that differ- 
ent components of Ta are naturally related to conditions 
on the lapse and shift respectively: 

EE fTa = N-^{dtN - N'd^N) + K, (46) 

+^'^r,,,g^\ (47) 

where '■^-'ryfe is the Christoffel symbol associated with 
the three-metric gij . We see that depends on the time 
derivative of the lapse, and that Ti depends on the time 
derivative of the shift. It is natural then to impose slic- 
ing conditions using the F^ component of the target gauge 
source function, and to impose shift conditions through 
the spatial components Fi. Once Fj. and Fi are speci- 
fied, the time component Ft is obtained from the identity 
Ft = NFf + N^Fk- Finally, we will want to express the 
target gauge source function in terms of the first-order 
GH Einstein system variables {tpab^Tiab, ^iab}, therefore 
the expression for Fq from Eq. (|19p will be useful: 

Fa = 5'-'*ya + t^Hfea - - UalP^'^Ii^,. (48) 

The remainder of this section presents a list of target 
gauge source functions, fa, that describe commonly used 
gauge conditions in numerical relativity. Slicing condi- 
tions are described in Sec. IIII Al and shift conditions are 
given in Sec. IIII Bl 

A. Slicing Conditions 

One of oldest gauge conditions used in numerical rel- 
ativity is maximal slicing 6], where the trace of the ex- 
trinsic curvature of the t = constant hypersurfaces van- 
ishes: K = 0. More generally constant curvature slicings 
are sometimes used, K = Kq, where Kq is constant on 
each slice (but may be a specified function of time). This 
gauge condition can be written in the form Gj: = 0, where 

Gi = Ko-K = Ko- i.9'm,, - g'H-<i>,j,. (49) 

Using Eq. (gS]) with Eq. (gH) and (gH), we obtain 

Ft - -^t"t'Hafc-piXo + ^(pi-l)5'-'ny 

-H(pi-l)tV^$.ja. (50) 

The choice of the arbitrary slicing gauge parameter, pi = 
1, gives a very simple expression for the constant curva- 
ture target gauge source function F^, but other choices 
may be more stable or more effective. 
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Perhaps the most widely used shcing conditions are 
various members of the family introduced by Bona and 
Masso 0- These gauge conditions are evolution equa- 
tions for the lapse N having the general form^ 



dtN^-N^f{N){K-Ko), 



(51) 



where f{N) is an arbitrary function of the lapse. The 
particular case f{N) — 2/N corresponds to the widely 
used one-plus-log shcing conditions [1, H, [13, UH- An 
expression for the general form of this gauge condition in 
terms of the first-order GH fields is given by 



V{N) 



1 



2Nf{N) 



NH''t''^,ab- (52) 



Using this condition in Eq. (|45|) results in the needed tar- 
get gauge source functions for these Bona-Masso slicing 
conditions: 

The slicing gauge parameter choice pi = 1 makes this ex- 
pression for the Bona-Masso gauge condition particularly 
simple; however, any choice with pi 7^ is allowed. 



B. Shift Conditions 

The simplest shift condition (from our perspective) is 
referred to as F- freezing This condition fixes the 

trace of the Christoffel symbol associated with the con- 
formal spatial metric gij = g^gtj , where g = det gij and 
A is a constant. (Often A is chosen to be A = — ^ so that 
det gij — 1, but any value is allowed.) The relevant trace 
of this conformal connection is defined by 

(54) 

The F-freezing shift condition simply requires that 

dt^^^r = 0. (55) 

For our purposes this must be translated into a condi- 
tion on the spacetime metric and its first derivatives. 
This is accomplished by integrating Eq. (j55p to obtain 
(3)fi ^ (3)fi(o)^ where (3)f*(0) is the trace evaluated at 



^ The original gauge condition in Q contains a derivative along 
the timelike normal instead of a partial time derivative. 



the initial time. This condition can be expressed in terms 
of first-order GH fields as 



' <f,M. (56) 



Using Eq. ([i5|) and this gauge condition is easily 

transformed into the needed target gauge source func- 
tion: 



- i[l - P2(l + A)].g^'=$,yfc - ^t''t''^,ab - t"n„ 



-P2g 



^g.,^'^t^0) + {p2-l)g'''<i>,M, 



(57) 



where the shift gauge parameter p2 ^ can be chosen 
freely. As a modest generalization we might also want to 
consider F-fixing conditions for which ^^'>T'^ is specified 
as a function of time. For example we might want to set 
(3)f^(i) = (3)f '(O)e-'''. This can be done by replacing 
(3)f *(0) with the desired (3)f *(t) in Eqs. ^ and (jFfl) . 

The most commonly used shift conditions in the nu- 
merical relativity community are the F-driver conditions. 
The simplest of these can be written as the following evo- 
lution equations for the shift (loj . 



dtN' = B\ 



(58) 
(59) 



where (3)p* jg the trace of the conformal spatial connec- 
tion, Eq. ([M]) . and v and 772 are adjustable constants. 
The parameter v is usually set to 1/ = | on the basis 
of causality arguments d, ■ But these arguments do 
not apply when the lapse and shift are evolved with the 
GH Einstein equations, so we leave v as an adjustable 
parameter. Unfortunately this shift condition is not of 
the form Gi — Gi{x,ip,dip), which is required by our 
gauge driver system, because the right side of Eq. ([59]) 
depends on second derivatives of the spacetime metric. 
This particular F-driver condition, Eqs. ([55]) and (|59p . 
can be transformed however into the more useful form 



(60) 
(61) 



We note that Eqs. ([60|) and ((6T|) are equivalent to 
Eqs. ([58|l and ((59|) when 772 7^ 0. This can be seen by 
differentiating i?* = z^f'^^F* — 772 T*] with respect to time 
to determine that Eq. ([55)) is equivalent to Eq. (|6ip. 

This transformed F-driver condition does not depend 
on the second derivatives of the spacetime metric, so it 
is of the form required for our gauge driver system. This 
F-driver condition, Eq. (|60p . can be written in terms of 
the first-order GH fields as 



N 



N2g^ 



1 + A 



9t3 



9i 9- 



5.V')$,fch (62) 
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where the auxihary field T' must be evolved using 
Eq. (pT|) . and is treated as an independent dynamical 
field along with the GH and gauge driver fields. The 
addition of Eq. (pT|) to the evolution system does not 
affect hyperbolicity. (The combined system has the ad- 
ditional characteristic fields T% all of which have charac- 
teristic speed zero.) When evolving Eqs. (|58p and 
it is common practice to set dtN^ = initially J^; the 
equivalent condition in our notation is initially choosing 
7^2 T' = '^^^r*. The target gauge source function Fi for 
this F-driver is obtained from Eq. (p^ using Eqs. 
and (gg): 



iV2 



(63) 



The shift gauge parameter p2 can be chosen to have any 
non-zero value. 



IV. FLAT SPACE STABILITY ANALYSIS 

The analysis of the gauge driver equations in the previ- 
ous sections is concerned with rather general questions, 
such as: Are the equations hyperbolic? What are the 
appropriate boundary conditions? How are particular 
gauge conditions implemented? In this section (and the 
next) we focus on questions about the stability and effec- 
tiveness of the gauge driver equations, such as: Are the 
gauge driver evolution equations stable? How well do 
the equations actually drive Ha toward the target gauge 
source function Fa7 In this section we use (mostly) an- 
alytical methods to explore these questions for simple 
cases that can be described as linear perturbations of flat 
spacetime. We consider three successively more compli- 
cated versions of this flat spacetime problem: First, we 
analyze the solutions to the gauge driver equation with 
a flxed Fa on a flat background spacetime. Second, we 
generalize this problem by allowing Fa to have a pre- 
scribed time dependence. Third, we analyze the more 
realistic case of the coupled gauge driver and GH Ein- 
stein systems for linear perturbations of flat spacetime. 
We present this analysis in some detail for the case of a 
target Fa representing Bona-Masso slicing and a F-driver 
shift condition. 

Before we specialize to these three specific problems 
however, we first establish some common notation and 
present the basic equations. Since we are perturb- 
ing about flat spacetime, it is convenient to decom- 
pose the solutions into spatial Fourier basis functions. 
Thus we assume that the spatial dependence of each 
of the perturbed fields is e'^''^^\ where kj is a (con- 
stant) wave vector, and are the spatial Cartesian 
coordinates. We assume that the gauge source func- 



tion Ha, the target function Fa and the time averag- 
ing field 9a have the forms Ha{t,x) = 6Ha{t)e^^'^\ 
Fa{t,x) = 5Fa{t)e'''^'=\ and 0a{t,x) = 59a{t)e'''^^\ We 
also assume that the spacetime metric ^at has the form 
tpab{t,x) = rjab + S^ab{t)e^'^''^\ whcrc rjab is the fixed 
background Minkowski metric with N = I, gij = Sij 
and A^* constant. With these assumptions the linearized 
gauge driver system, Eqs. ([^ - ([^ and (P7)l . can be 
written in the form: 



dtSHa 



dtSU 



H 



i(3kSHa = 

-SU^ -J^N^{S^fa-tk,SHa), 

[2^2(1 - 6) - «/3fc]<5Hf = -ik^<P^a 

+m5ea+p.lil-^l)iSHa~6Fa) 

-^f-i^W{5^fa-^k,5Ha), 

dt{S<ffa - tkjSHa) ^ 

+ikj-f^N\5^^a-^h5Ha) 
-{^^ -ipk){5^fa-^k,5Ha), 
mS0a = [2p2tkf3{l - Cs) + P-lii]6Ha 
+ [2m2(6 - 6) - iPk] + ik'S^fa 
-2/12(1 - ^sh^N^iS'^l - ^kjSHa) 



(64) 



(65) 



(66) 



ikjSHa 



^lUl5Fa, (67) 



where k^ = k^kj, 13k = kjN^ , and contractions are done 
with the fiat background metric. 

This Hnearized gauge driver system, Eqs. ([M1) ~ ([S7)) . 
can be simplified somewhat. We note that Eq. ([55]) 
implies that violations in the gauge constraint SC^a = 
5^fa ~ ikjSHa always decrease toward zero exponentially 
on the timescale 1/7|^. Since the system is linear, we 
can (without loss of generality) confine our attention to 
the constraint satisfying solutions, 5^fa — ikjSHa- This 
condition and Eq. can be used to eliminate the fields 
SU-a and S^fa from the system, resulting in the following 
simplified evolution system for 5Ha and 59a'- 

d^SHa + 2 [/X2 (1 - 6) - if3k] dtSHa 

-f [fc2(l - /?2) + fj^l^i _ _ 2z/i2/3A:(l - 6)] SHa 
^~m69a + P-lil-^i)SFa, (68) 

dt69a + mSda = - [2/i2(6 - 6) - iPk]dt6Ha 



fc2(l-/32) 
^^^^lSFa. 



2i^^2|3k{l - 6) 



SHa 



(69) 



We note that the linearized gauge driver system, 
Eqs. (|68|) - ((69)) . does not depend on the metric pertur- 
bations Sipab, except through the target function SFa- 
This rather weak coupling means that the gauge driver 
equations just respond to whatever target SFa the gauge 
and spacetime geometry dictate. It makes sense then to 
investigate the intrinsic response of the gauge driver sys- 
tem to a given SFa- We consider two simple test cases. 
First, in Sec. IIV Al we consider the case where the target 
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gauge source function, 5Fa, is time independent. Second 
in Sec. lIVBi we consider the more general case where 
6Fa = SFa{t) is a prescribed function of time. Then fi- 
nally, in Sec. IIV CI we consider the more interesting and 
realistic case where the gauge driver and GH Einstein 
systems are coupled, using the target SFa appropriate 
for Bona-Masso slicing and a F-driver shift condition. 



A. Time Independent SFa. 

We consider first the case where the target gauge 
source function has the form. Fa = SFae^''^^\ for con- 
stant SFa- We also assume that the shift of the back- 
ground spacetime vanishes: (3 = 0. In this case the gen- 
eral solution to Eqs. ((68|l and ([69| has the form: 



6Hait)=5Fa + J2SH:e'"\ 



%{t) = --SFa 



E 



P + 2s^/i2(g3-6)-Mi6 

Sr + -111 



(70) 



SH:e^^\ (71) 



where the SHa are constants and Sr are the roots (as- 
sumed to be non-degenerate) of the characteristic poly- 
nomial, 

o = 4 + [2fi2{i~C2) + m]4 

+ [e + - Ci) + 2/i2f?i(l - 6)1 Sr + mM?- (72) 



Equation (|72p is the necessary and sufficient condition 
that the solution satisfies Eqs. ([55)1 and ([55)1 . The three 
roots of Eq. ((72|) consist of a real root, so, and a complex 
conjugate pair, s±. 

Figure [1] illustrates the dependence of the real parts of 
the roots, so and s±, on the wavenumber k for the case 
/i = /ii = /i2 = '71 and = = ^2 = Cs- These roots 
have strictly negative real parts for all k, so the gauge 
source function 5Ha is always driven toward the target 
gauge source function 5Fa- At least for this simple case, 
SHa approaches the target SFa exponentially. 

Simple analytical expressions for the roots of the char- 
acteristic polynomial, Eq. (|72p. exist in the limits of small 
and large k. The large k limit is the most interesting, 
because it describes the sufficently short wavelength per- 
turbations of any spacetime. The asymptotic expressions 
for the large k roots are, 



Re (so) 
Re(s±) 



(73) 



(74) 



These results show that the s± modes are damped at 
approximately the rate ^2(1 — C2) + ?yi/2 in the large k 
limit, while the damping rate for the so mode approaches 



0.0- 



-0.5 



1.0 



■1.5 



1 1 1 1 1 1 1 1 1 


- / 


ReCs^iLi) : 




Re(s^/|i) : 


1 . 1 . 1 . 1 . 



2 4 6 8 10 
k/|Li 



FIG. 1: Real part of the characteristic frequencies of the gauge 
driver system: so and s±. 



zero. These modes are stable for large enough fc, then, 
as long as rii^\ > and ^2(1 — C2) + ??i/2 > 0. 



B. Time Dependent 6Fa 



Next we consider solutions to Eqs. (p5|) and for 
the case where SFa is a specified function of time: SFa = 
SFa{t). In principle the equations could be solved ana- 
lytically by Laplace transforming the equations in time, 
and solving for each frequency component of 5Ha{t) sep- 
arately. Instead it is more straightforward, and perhaps 
more instructive, to integrate the equations numerically 
for some ihustrative SFait). We assume for this simple 
example that the shift of the background spacetime van- 
ishes, /3 = 0, and the other parameters that determine the 
system take the values: /i = /ii = /i2 = ?7i and = ^1 = 
^2 = ^3- We have solved the resulting simplified equa- 
tions numerically for the case SFa{t) = 3 + e^^*^^°^ 
with fc = 1. Equations ([55]) and ([55]) require initial con- 
ditions for SHa, dtSHa and S9a. We use SHaiO) = (5FJ0), 
dtSHaiO) = and fxS0a{O) = -k'^SHa{0). These initial 
data for SHa and its time derivative were chosen to be 
fairly well matched with the target SFa- They are similar 
to the initial conditions used in our more realistic tests 
in Sec.|Vl This target SFa changes significantly for times 
near t = 10, so this test explores how well the gauge 
driver system is able to track an evolving target SFa- 
Figure [2] shows that the gauge driver equation is fairly 
successful (at the few percent accuracy level) in driving 
SHa{t) toward SFa{t) for /i > 2 even in this rather dy- 
namical situation. 
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FIG. 2: Response of the gauge driver system to a time depen- 
dent SFa of the form, SFa = S + e"'*"^"' initial conditions 
SHa{o) = SFa{0), dt5Ha{0) = 0, and a range of values for the 
damping parameter /i = /ii = /i2 = This test uses k = 1 
and = 6 = = 0. 



C. Coupled Systems 

Finally we investigate the stability of the coupled 
gauge driver and GH Einstein equations for perturba- 
tions of flat spacetime. The perturbed Einstein system 
reduces to a relatively simple form"^ in this case: 



r^'^dcdMab + da5Hb + dbSHa = 0. 



(75) 



We study the stability of the coupled system, Eqs. 
(j69p , and (|75p , by Laplace transforming the equations in 
time, i.e., by considering solutions with time dependence 
e**. In this case Eqs. ([55]) and can be reduced to the 
single equation, 



where P{s) is defined by. 



5 Fa. 



(76) 



p{s) - f + 2^i2{l-^2)s + k'' + 
m 



\e - - 2^/3A;^2(l - 6) 
+s[2M2(6-6)-«/?fc]}- (77) 



^ In this analysis we assume that the gauge constraint SHa + S^a = 
is satisfied. The analysis of the GH Einstein constraint evo- 
lution system in Ref. Q] shows that violations of this constraint 
are damped exponentially for perturbations of flat spacetime. 



We use the notation s = s — ikf3. The analogous ex- 
pressions for the Laplace transform of the GH Einstein 
system, Eq. (|75p . are given by 







= (J^ + A:^)^Vt 



sSHj — ikjSHf, 



= {s^ + k^)Sijji - ikjSHi - ikiSHj, 



(78) 
(79) 
(80) 



where SHj: — SHt — N^SHj, etc. These equations can be 
used to express 6Ha and Sipji in terms of S'tpfa the 
case s 7^ 0: 



5H, 



2s 



s 

s-2 



(81) 
(82) 



s ^ (^iskjSipjii + iskiStpfj + kjkiS^p^^) . (83) 



The case s 

U2 



case is essentially trivial: In this case 
k'^Sipff — 0, fc^Vtj = ikjSHf, k'^Sipji — ikjSHi + ikiSHj 
and 5Ha = SFa- The metric perturbation in this case 
is pure gauge (an infinitesimal coordinate transforma- 
tion generated by the time independent SHa/k'^), and 
the gauge source function SHa is identical to the target 
SFa in this case. So we focus on the s ^ case for the 
remainder of this discussion. 

We consider in detail now the coupled gauge driver 
system for the case of Bona-Masso slicing and the T- 
driver shift condition. The perturbed flat space limit of 
SFa for the Bona-Masso driver, Eq. (|53p. is given by 



SFi 



/(I) -pi ipkpi 



2/(1) 2/(1) 

+t{pi - l)k^St, 



Hit 



l)sS^b'i, (84) 



while the target for the T-driver shift condition, Eq. (j63p . 
reduces to 



SFj = (s - p2s)Stp^ 



VP2S 

s + m 



1 Si,jik' 



-—kiS^ff 

2 3 1 tt 2 



vp2s(l + A) 



1 



kjS^\. (85) 



The spatial metric perturbations, Sijjji, that appear in 
Eqs. dll]) and can be replaced by S^fa using Eq. ((83|l : 



SFi 



/(l)-pi i(3kpi P{pi~l) 



2/(1) 



SF, = - 
^ 2 



fc2 

— I J^P2S 

V s 



2/(1) 
1 - A 



2s 



Vp2S 



- 1 



-1-1 



s \s H 

VP2SX 



kjk^Sipii- 



P2S - s 



Hu 



Now substitute these expressions for SFa , Eqs 
(I57)l . and the expressions for SHa, Eqs. ([5T|) and 



(87) 
and 
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into the perturbed gauge driver Eq. ([76]) . The result 
is a system of hnear algebraic equations for Stp^^. This 
system can be decoupled and non-trivial solutions exist 
if and only if the frequency s satisfies one of the following 
characteristic polynomials: 







P{s)+^^i{ 1 



s + l]l 



Pi - /(I) , ifikpi , , P 
s — h + [pi - 1) — 



/(I) 

'-P{s)+^ii\ 1 



/(I) 



S + 7]l 



kW 1 - A 
— iyp2s— 



P{s) 



VP2S 

s + r]2 



- 1 



P2S - 



P2S - S 



(88) 



(89) 



(90) 



where P{s) is defined in Eq. (|77p. 

The flat space stability analysis presented here is rele- 
vant to generic spacetimes when the wavenumber k of the 
perturbation becomes sufficiently large. We have solved 
the characteristic polynomials in Eqs. ([88|) -([90 |) in this 
limit. The leading order expressions for the real parts 
of these roots are given as follows. For the time slicing 
modes (in which (^-044 0)? the roots of Eq. ([88|) . we have 



Re(s) 
Re(s) 



mPlPi 



o{k-% 



(l-/32)2/c2 

±i| [771 + 2^2(1 -6)] 



(91) 



-4Ai?pi(l-Ci) 



l"/(l)±/3 
/(I) 



1/2 



[771+2/12(1-6)] +0 (A:-') 



(92) 



The asymptotic forms of the roots of the longitudinal 
modes (in which k^S^^j ^ 0), Eq. ([89|) . are given by 



Re(s) = -f]2 + 0{k-^), 

Re(s) ±i{[r,i + 2M2(l-6)]' 

-i [,71 + 2/^2(1-6)] +Oik-'). 



(93) 



1/2 



(94) 



Finally the asymptotic forms of the roots of the trans- 
verse modes (in which [/c^g*-' — k^k^]6^/jfj ^ 0), Eq. ([gO)) . 
are given by 



Re(s) = -rj2 + 0{k-^), 

Re(s) = ±i{ [771 + 2/^2(1 -6)]' 



(95) 
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FIG. 3: Maximum damping rate of the modes as a function of 
the Bona-Masso slicing condition parameter /(I). The other 
parameters used for this case are k = 1, l3 — 0, p, = pi = 
M2 = 771 = ^7?2, 6 = 6 = 6 = 0, pi = P2 = = |, and 



-4M?P2(i-ei)[i±/3-j^]} 



1/2 



--[771 + 2/72(1-6)] +0(fc"')- (96) 

All four ± sign combinations represent distinct roots in 
Eqs. (HH), dM]), and ([M]). Stability of the gauge driver 
system requires Re(s) < 0. Therefore, stability of the 
short wavelength modes requires the following inequali- 
ties on the system parameters: 



< 771/71, 

< 771 + 2/^2(1-6), 

l-/(l)±/9 



< Pi(l-6)- 



/(l) 



< ?72, 

< P2(l-6)[l±/5-i^(l-A)], 

< P2(i-6)[i±/3-^^]- 



(97) 
(98) 

(99) 

(100) 
(101) 
(102) 



We note thate these conditions can be satisfied for small 
values of f3 by taking 7/1 > 0, 772 > 0, p2 > 0, 6 < 1, 
6 < 1, pi > 0, P2 > 0, < /(I) < 1, < 1, iy{l-X) < 1. 

We have also explored the roots of these characteristic 
polynomials numerically. Figure[3]illustrates Max[Re(s)], 
the root of these equations having the largest real part, 
as a function of the parameter / (I) that characterizes the 
Bona-Masso slicing condition in this fiat space limit. The 
curves correspond to the roots for the driver system with 
various values of p = pi = P2 = Vi = gi'?^, = 
1/ = |, and k = 1. These parameter values were chosen 
because they satisfy the inequalities in Eqs. (|571) - (|102p . 
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FIG. 4: Maximum damping rate of the modes as a function of 
the wavenumber k. The other parameters used for this case 
are /? = 0, = /^i = ^2 = ?7i = ^V2, ^1=^2 = ^3= 0, 
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FIG. 5: Maximum damping rate of the modes as a function 
of the shift parameter (3. The other parameters used for this 
case are = 1, ^ = ^1 = ^2 = J?! = ^')2, = ^2 = ^3 = 0, 
Pi = P2 = |, /(I) = i, = |, and A = 



and because they perform fairly well for the 3D numerical 
tests discussed in Sec. |Vl We see that the maximum real 
part of s is negative for /(I) in the range < /(I) < 
1, and so the coupled gauge driver system is stable for 
these values. The system is most stable for /(I) w i, so 
we adopt this value in our numerical tests of the gauge 
driver system in Sec. |Vl We also note that the standard 
value, /(I) = 2, used for one-plus-log slicin g b y most 
of the numerical relativity community [1, lol. lid. HH is 
unstable when used in our gauge driver equations. This 
does not imply that /(I) = 2 is a bad choice when used 
in a standard three-plus-one evolution, only that it is 
unstable when used with our gauge driver system. 

Figure S] illustrates the k dependence of Max[Re(s)] 
for a range of values of the gauge driver damping coeffi- 
cients fJ- ^ 1^1 = IJ.2 = Vi — ^^2- For short wavelength 
perturbations, i.e. for values of k with k > ji, Max[Re(s)] 
decreases as fi increases. Thus the solutions with large 
k are damped more effectively as ^ increases. However, 
for long wavelength perturbations, i.e. for values of k 
with k < fj,, Max[Re(s)] increases as /i increases. Thus 
the solutions with small k are less efficiently damped as 
IJ, increases. It follows that there is an optimal value of /i 
to use for any particular problem: choose /x « kc, where 
1/kc corresponds to the lengthscale on which the gauge 
condition needs to be enforced most effectively. 

Figure [5] illustrates the dependence of Max[Re(s)] on 
the background shift parameter (3 for a range of values 
of the gauge driver damping coefficients fj, = ni = fj,2 — 
rji — ^772 . For small values of (3 we see that the system is 
stable, however for /3 > ^ the system becomes unstable. 
This instability may be important in more realistic prob- 



lems that involve black holes. Even for single black-hole 
spacetimes, the usual time-independent coordinate rep- 
resentations have non-vanishing shifts with /3 « 1 near 
the horizon. Binary black hole spacetimes also use large 
shifts (with /3 > 1 in many cases) when coordinates that 
co-rotate with the black holes are used. We explore the 
stability of this gauge-driver system for the case of single 
black-hole spacetimes in Sec. |Vl 

We have also examined several other slicing and shift 
conditions using these perturbation techniques. As a 
consequence of section IIV Al our gauge driver system is 
stable for harmonic slicing 6F^ = and harmonic shift 
SFi = 0. We also find that the combinations of a stable 
Bona-Masso slicing condition with harmonic shift, and 
of harmonic slicing with the F-driver shift condition are 
stable. However, we find that the maximal slicing and 
F-freezing conditions are unconditionally unstable when 
enforced through our gauge driver equations. 



V. NUMERICAL TESTS 

In this section we describe the results of 3D numerical 
tests of the gauge driver system. We consider two cases: 
first a Schwarzschild black hole with perturbed lapse and 
shift, and second a Schwarzschild black hole with a super- 
imposed outgoing physical gravitational wave pulse. The 
full coupled non-linear GH Einstein and gauge driver sys- 
tems are solved numerically for these cases. We measure 
the stability and effectiveness of the gauge driver system 
in these tests as it attempts to drive the gauge toward 
Bona-Masso slicing and F-driver shift conditions. 
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These numerical tests are conducted using the infra- 
structure of the Caltech/Cornell Spectral Einstein Code 
(SpEC). This code uses pseudospectral collocation meth- 
ods, as described for example in Refs. fl^, [l3|. We use 
the generalized harmonic form of the Einstein equations, 
as described in Ref. The evolution equations for the 
combined GH Einstein and the gauge driver system are 
integrated in time using the adaptive fifth-order Cash- 
Karp method [lH| ■ We use a form of spectral filtering, as 
described in Ref. [l3] , that sets to zero in each time step 
the changes in the top four tensor spherical harmonic 
expansion coefficients of each of our evolved quantities. 
This filtering step is needed to eliminate an instability as- 
sociated with the inconsistent mixing of tensor spherical 
harmonics in our approach. 

Initial conditions are needed for any evolution of the 
combined GH Einstein and gauge driver systems, and 
these initial data consist of the spacetime metric ipab, 
the gauge source function Ha, and their time derivatives. 
For the tests described here we take the initial spacetime 
metric ipab to be the Schwarzschild geometry plus per- 
turbations as described in Secs. lV Al and fVBl We set the 
time derivatives of the spatial components of the metric 
to zero for the pure gauge perturbation test in Sec. IV Al 
and equal to the appropriate time derivative of the su- 
perimposed physical gravitational wave pulse for the test 
in Sec. IV Bl The remaining initial data needed for these 
evolutions, dtN, dfN^, Ha, and dfHa, are pure gauge 
quantites. The time derivatives of the lapse and shift are 
chosen here to ensure that Ha satisfies the desired gauge 
condition. Ha = Fa , initially. And finally the initial value 
of Ha is chosen here to ensure that the gauge constraint, 
Ca ~ Ha + Ta = 0, vanishcs initially. 



A. Black Hole with Gauge Perturbation 

For this test we consider a Schwarzschild black hole 
with perturbations in the lapse and shift. For the un- 
perturbed hole we use isotropic spatial coordinates and 
maximal time slices (l6l . [l7| . The unperturbed spatial 
metric in this representation is given by. 



where 



g.jdx'dx^ = ^-^ (dx^ -I- dy'^ + dz^) , (103) 



y 



, and R{r) (the areal radius) 



satisfies the differential equation, 



dR _ R I 2M 
dr ~ 7V ~ ~R~ 



C2 
R^' 



(104) 



The constant M is the mass of the hole, and C is a param- 
eter that specifies the particular maximal slicing. Finally, 
the unperturbed lapse N and shift N'^ for this represen- 
tation of Schwarzschild are given by. 
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(106) 



where f' is the outward directed radial unit vector: 
^ij-PfJ = 1. 

We perturb this spacetime by changing the initial val- 
ues of the lapse and shift, and their time derivatives. This 
type of perturbation changes the spacetime coordinates 
(or gauge) of the solution, but not its geometry. For 
these tests we perturb the lapse and shift of Eqs. (jlOSp 
and (|106|) by adding, 

5N = Asin(27rr/ro)e-('-'''=)'/™V/™, (107) 

5N' = ^sin(27rr/ro)e-('-'^<=)'/"'V/™f\ (108) 

where Yim is the standard scalar spherical harmonic. In 
our numerical tests we use the background metric with 
C = 1.73M2, and perturbations with A = 0.01, = 
15M, w = 3M, I = 2, m — 0, and various values of the 
radial wavelength tq. 

We perform these numerical tests on a computational 
domain consisting of a spherical shell that extends from 
r = 0.78Af (just inside the horizon in these coordinates) 
to r = 30M. We divide this domain into eight subdo- 
mains. In each subdomain we express each Cartesian 
component of each dynamical field as a sum of Cheby- 
chev polynomials of r (through order Nr — 1) multiplied 
by scalar spherical harmonics (through order L). The 
radii of the inner and outer edges of the various sub- 
domains are adjusted to distribute the truncation error 
for this problem more or less uniformly. The specific 
radii of the subdomain boundaries used in this test are 
0.78M, 2.38M, 4.6M, 8.83M, 13.07M, 17.30M, 21.53M, 
25.77M, and 30. OM respectively. The values of the pa- 
rameters associated with the gauge driver system used 
for this test are: v = |, A = — |, pi — p2 — \, 
^1 = ^2 = ^3 = 0, and various values of the parameter 
/i = = /i2 = ?7i = ■^'n2- The Bona-Masso slicing con- 
dition includes a target value for the extrinsic curvature 
Ko; for this test we set = 0. 

Figure [H] illustrates the constraint violations for a set 
of representative evolutions from this test, and demon- 
strates the exponential convergence of our numerical 
method. The solid curves represent the constraints as- 
sociated with the GH Einstein system, while the dotted 
curves represent the constraints of the gauge driver sys- 
tem. We measure the constraint violations of the GH 
Einstein system for these tests using the norm 1 1 Cqh | | de- 
fined in Eq. (71) of Ref. The norm 1 1 Cqh | | is scaled so 
that it becomes of order unity when constraint violations 
start to dominante the solution. We define an analogous 
norm 1 1 Ch 1 1 for the gauge driver system: 



|Ch| 



^m"'g'^ (d,H,d,Ha -f 9.Hf 9,Hf 



(105) 



.(109) 
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FIG. 6: Solid curves show the constraints of the GH Ein- 
stein system 1 1 Cgh | | , dotted curves show the constraints of 
the gauge driver system || ChI| for a test with — 1.0/M and 
radial wavelength ro = 6.0M evolved at several numerical 
resolutions. 

The quantity to"'' is a positive definite matrix, which we 
set to the identity matrix, m'^'' = S'^'', for these tests. Fig- 
ure [6] shows the constraints for a particular test run with 
/i = 1.0/M and ro = 6.0M. The analogous graphs for the 
other tests reported here are qualitatively similar, with 
somewhat larger but still convergent "spikes" in 1 1 Cqh \ \ 
at early times (t < 25M) for the ^ = 0.5/M case. Fig- 
ure [H] shows that the constraints are well satisfied in our 
evolutions, and demonstrates that our numerical meth- 
ods are (exponentially) convergent. The mild power law 
growth in the constraints seen at late times is sublinear, 
and is not something that concerns us. 

Figure[7]illustrates the effectiveness of the gauge driver 
system, at least for this test problem. We measure the 
difference between the gauge source function Ha and the 
target function to which it is being driven, Fa, using the 
following norm: 

\\H-F\\^ ^ J^m-\Ha-Fa)iH,-F,)d^x 

where (as before) the matrix to"'' is set to the identity, 
^ab _ jab^ thcsc tcsts. This norm is scaled so that 
Ha bears little resemblance to the target Fa whenever the 
norm becomes of order unity. Figure [7| shows that the 
gauge perturbation used in this test violates the desired 
gauge conditions rather severly at early times. The norm 
I |i? — i^l I/I |F| I is driven to values as large as 0.7 at about 
t = 20M when the ingoing part of the gauge perturbation 
interacts most strongly with the black hole. After this 
initial interaction, the gauge driver system takes over and 
effectively drives \ \H — F\\/\\F\ \ to values below 10"'^ on 



FIG. 7: Effectiveness of the gauge driver equation is demon- 
strated by showing \\H — F\\/\\F\\ for evolutions with ra- 
dial wavelength ro = 6M and several values of the gauge 
damping parameter fiM £ {i, 1, |,2}. These tests evolve a 
Schwarzschild black hole with strongly perturbed lapse and 
shift. 



timescales of 40Af to 60A/, depending on the value of the 
gauge damping parameter /i used in the evolution. Fig- 
ure[7]shows evolutions of gauge perturbations with radial 
wavelength rg = 6A/, and several values of the damping 
parameter jiM G {i,l,|,2} computed with numerical 
resolution A^,, = 17 and L = 13. The gauge driver sys- 
tem is more effective at reducing ||i7 — F||/||i^|| quickly 
at early times, t < 75M, for larger values of /i. 

The fi = 2/M case shown in Fig. [7] has a mild in- 
stability, that first appears at about t — 300M. This 
is a gauge instability since it does not affect any of the 
constraint quantities. Larger values of /i are progres- 
sively more unstable. This instability may be related 
to the rather unusual dispersion relation for this gauge 
driver, as shown for the flat space case in Fig. [D The 
gauge driver equation becomes increasingly ineffective 
for driving the long wavelength components of Ha to- 
ward Fa as /i increases. This poor damping efficiency 
for long wavelengths, together with our rather simplistic 
boundary conditions or the inherent instability associ- 
ated with large shifts, may well be the cause of this in- 
stability. Figure [5] provides some additional insight into 
the way the gauge driver equation responds to differ- 
ent perturbations. The evolutions shown in Fig. [8] are 
all performed with fi — 0.5/M but several different val- 
ues of the radial wavelength of the gauge perturbation: 
ro e {4M,6M,8M, lOM}. These tests show that the 
gauge driver system causes ||7f — i^||/||F|| to approach 
zero more quickly (at least at early times) for shorter 
wavelength perturbations. Our hope is that this ability to 
efficiently control short wavelength features of the gauge 
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FIG. 8: Effectiveness of tfie gauge driver equation is demon- 
strated by sliowing \\H — for evolutions witli n = 
0.5/M and several values of the radial wavelength of the per- 
turbation ro e {4M, 6M, 8M, lOM}. These tests evolve a 
Schwarzschild black hole with strongly perturbed lapse and 
shift. 



FIG. 9: Effectiveness of the gauge driver equation is demon- 
strated by showing — for evolutions with /i — 
0.25/M obtained with a variety of numerical resolutions. This 
test uses a Schwarzschild black hole with a superimposed out- 
going gravitational wave pulse. 



is what will be needed to prevent the kinds of localized 
gauge singularities that often appear in our evolutions of 
binary black hole spacetimes. 

B. Black Hole with Physical Perturbation 

Our second numerical test of the gauge driver system 
uses a Schwarzschild black hole with a superimposed out- 
going gravitational wave pulse, as described in Refs. [H, 
|l8|. The background solution is a Schwarzschild black 
hole in Kerr-Schild coordinates, 

2M 

ds^ = ~dt^ H {dt + drf + dx^ + dy^ + dz^, (111) 

r 

where r"^ = x"^ + y'^ + and M is the mass. We super- 
impose an odd-parity outgoing quadrupolar gravitational 
wave perturbation constructed using Teukolsky's method 
p^ . Its generating function is taken to be a Gaussian, 
G(r) = ^exp[-(r - rc)^/w;^], with ^ = 4 x 10"^, Tc = 
5M, and w = 1.5M. Using this perturbed Schwarzschild 
solution as the input conformal metric, the full non-linear 
initial value equations (in the conformal thin sandwich 
formulation) are solved to obtain initial data that sat- 
isfy the constraints [l^]. This procedure yields initial 
values for the spatial metric, extrinsic curvature, lapse, 
and shift. We note that the resulting solution to the 
constraints is still nearly (but not completely) outgoing. 

The computational domain for this test problem is 
taken to be a spherical shell extending from r = 1.9Af 
(just inside the horizon in these coordinates) out to 



r = 41.9A{f. This domain is subdivided into four 
spherical-shell subdomains of width Ar — lOAf . On 
each subdomain, the numerical solution is expanded in 
Chebyshev polynomials and spherical harmonics as be- 
fore. For these tests we use numerical resolutions with 
€ {21,31,41,51} coefficients per subdomain for the 
Chebyshev series and / ^ L with L e {8, 10, 12, 14} for 
the spherical harmonics. 

Figure ini illustrates the effectiveness of the gauge driver 
equation for imposing the Bona-Masso slicing and F- 
driver shift conditions in evolutions of a Schwarzschild 
black hole with physical gravitational wave perturbation. 
These tests were performed with the gauge damping pa- 
rameter /i = 0.25/M. For this test we set the target 
value for the extrinsic curvature Ki^ to that of an un- 
perturbed Kerr-Schild spacetime. The various curves in 
Fig.iniillustrate how _F||/||F|| changes for evolutions 
performed with different numerical resolutions. The re- 
sults are qualitatively similar to those of the first test: 
the black hole with physical gravitational wave perturba- 
tion does not satisfy the target gauge conditions exactly 
at early times, but the gauge driver equation reduces 
— i^||/||F|| to very small values by about t = 75 Af. 
This test is less severe in some sense than our first pure 
gauge perturbation test, since the initial data in this case 
contains an outgoing gravitational wave pulse that never 
interacts very strongly with the black hole. 
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VI. DISCUSSION 

We have presented a new gauge driver evolution sys- 
tem in Sec. |lT] that makes it possible to impose a wide 
range of gauge conditions in the generalized harmonic 
(GH) formulation of the Einstein equations, without de- 
stroying its hyperbolicity. The key idea is to construct 
an auxiliary hyperbolic evolution equation for the gauge 
source function Ha that drives it toward the desired tar- 
get Fa . Section IIIII shows how many of the gauge condi- 
tions widely used by the numerical relativity community 
can be included in this way. In Sec. IIVI we analyze the 
effectiveness and stability of the combined GH Einstein 
and gauge driver system for the case of perturbations 
of flat spacetime. This analysis shows that the gauge 
driver equation effectively drives Ha toward F^, when Fa 
is specified a priori as a function of the spacetime coor- 
dinates. We were somewhat surprised to find, however, 
that the gauge driver system can be quite unstable when 
it is coupled to the GH Einstein system. We found that 
common gauge conditions like maximal slicing and the 
F-freezing gauge conditions are unconditionally unstable 
when implemented using our gauge driver equation. This 
does not imply of course that those conditions arc unsuit- 
able for use with other forms of the Einstein system (like 
BSSN), just that they cannot be implemented in a com- 
pletely stable way in the GH Einstein system coupled 
to the particular gauge driver equations introduced here. 
Fortunately, we were able to find some of the commonly 
used gauge conditions that can be implemented in this 
way: certain Bona-Masso slicing conditions and a com- 
monly used form of the F-driver shift conditions. Our 
3D numerical tests in Sec. IVl show that the gauge driver 
system can impose these gauge conditions stably and ef- 
fectively for the evolutions of perturbed single black hole 
spacetimes. 

There has been a great deal of discussion in the litera- 
ture about the formation of shocks when certain dynam- 
ical gauge conditions are imposed [U, [H, . However, 
these discussions do not apply when those same gauge 
conditions are imposed using a driver condition. The 
gauge driver system imposes the desired gauge condition 
only approximately, not exactly. At best, the desired 
gauge condition is imposed exactly only asymptotically 
in time as the system approaches a time independent 
equilibrium state, and even in this state shocks do not 
necessarily form. On the contrary, there are many so- 
lutions even to bad gauge conditions that do not have 
shocks. What determines whether an evolution system 
develops shocks is the structure of the operator that 
evolves the spacetime metric and auxiliary fields. Our 
evolution system (including the gauge driver system) has 



been carefully designed to be linearly degenerate, a con- 
dition that prevents the formation of shocks (resulting 
from a crossing of characteristics) from smooth initial 
data [2^. Linear degeneracy does not prevent the for- 
mation of curvature singularities, of course, or even the 
formation of coordinate singularities that may arise from 
non-linearities in the non-principal parts of the evolution 
equations. 

Causality is another issue that appears to be less re- 
strictive for our gauge driver system than it is for directly 
imposed gauge conditions. For example, the parameter ly 
that appears in the F-driver system discussed in Sec. lHIB] 
must take values in the range < < | in order for that 
F-driver to evolve the shift in a causal way in the BSSN 
system p5| . There is no such restriction on j/, however, 
when this F-driver is imposed through our gauge driver 
system. In our system the shift is evolved, along with the 
rest of the spacetime metric, by the GH Einstein system. 
This system is manifestly hyperbolic, and all of the fields 
propagate within the physical light cone, no matter what 
target gauge source function is used in the gauge driver 
system. 

It is easy to imagine that the system presented here 
could be improved in several ways. It may be possible, 
for example, to improve the performance of the system by 
formulating boundary conditions for Ha that impose the 
desired gauge condition Ha = Fa exactly at the bound- 
aries. It may also be possible to formulate a different evo- 
lution operator for Ha that drives it more stably and/or 
more efficiently toward the desired target Fa- Finally it 
may be possible to find better target gauge conditions Fa- 
The ones studied here are those which have been found 
useful in evolutions of traditional three-plus-one formu- 
lations of the Einstein system like BSSN. But there may 
exist gauge conditions having much better stability and 
effectiveness properties when used as target gauge condi- 
tions within a gauge driver system. These questions, and 
others, will be addressed in future work on this problem. 
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